main_df <- main_df %>% mutate(group.treatment = factor(group.treatment, levels = c("random", "merit", "patronage")))

any_vign_bivariate_1 <-lm(any_vign ~ factor(group.treatment), main_df %>% filter(pub_official_indicator == 1, !(session.code %in% c("q30x01sj", "qtj4xm8s"))))
all_vign_bivariate_1 <- lm(all_vign ~ factor(group.treatment), main_df %>% filter(pub_official_indicator == 1, !(session.code %in% c("q30x01sj", "qtj4xm8s"))))
avg_vign_bivariate_1 <- lm(avg_vign ~ factor(group.treatment), main_df %>% filter(pub_official_indicator == 1, !(session.code %in% c("q30x01sj", "qtj4xm8s"))))

any_vign_bivariate_2 <-lm(any_vign ~ factor(group.treatment), main_df %>% filter(pub_official_indicator == 1, (session.code %in% c("q30x01sj", "qtj4xm8s"))))
all_vign_bivariate_2 <- lm(all_vign ~ factor(group.treatment), main_df %>% filter(pub_official_indicator == 1, (session.code %in% c("q30x01sj", "qtj4xm8s"))))
avg_vign_bivariate_2 <- lm(avg_vign ~ factor(group.treatment), main_df %>% filter(pub_official_indicator == 1, (session.code %in% c("q30x01sj", "qtj4xm8s"))))

observations <- c(nobs(any_vign_bivariate_1),nobs(all_vign_bivariate_1),nobs(avg_vign_bivariate_1),
                  nobs(any_vign_bivariate_2),nobs(all_vign_bivariate_2),nobs(avg_vign_bivariate_2))

any_vign_bivariate_1 <- coeftest(any_vign_bivariate_1, vcov=vcovHC(any_vign_bivariate_1,type="HC0"))
all_vign_bivariate_1 <- coeftest(all_vign_bivariate_1, vcov=vcovHC(all_vign_bivariate_1,type="HC0"))
avg_vign_bivariate_1 <- coeftest(avg_vign_bivariate_1, vcov=vcovHC(avg_vign_bivariate_1,type="HC0"))
any_vign_bivariate_2 <- coeftest(any_vign_bivariate_2, vcov=vcovHC(any_vign_bivariate_2,type="HC0"))
all_vign_bivariate_2 <- coeftest(all_vign_bivariate_2, vcov=vcovHC(all_vign_bivariate_2,type="HC0"))
avg_vign_bivariate_2 <- coeftest(avg_vign_bivariate_2, vcov=vcovHC(avg_vign_bivariate_2,type="HC0"))

table <- list(any_vign_bivariate_1, all_vign_bivariate_1, avg_vign_bivariate_1, any_vign_bivariate_2, all_vign_bivariate_2, avg_vign_bivariate_2)

note_text <- paste("Beta coefficients from OLS regression. Standard errors were calculated using the Huber-White (HC0) correction. 
                   The outcomes measures are indices capturing whether public official offered Weberian responses to (1) any of the vignettes; 
                   (2) all of the vignettes; and (3) an average indexed measure of responses.")

table = stargazer(table, type = 'latex', 
                  title = "The Effect of Selection Mode on Weberianness",
                  label = 'tab:weber_session',
                  model.names = F,
                  model.numbers = T,
                  digits = 3,
                  column.separate = c(1,1 , 1, 1, 1, 1),
                  column.labels = c("Any", "All", "Avg", "Any", "All", "Avg"),
                  dep.var.labels = NULL, 
                  add.lines = list(c("Session", "1", "1", "1", "2", "2", "2"),
                                   c("Observations", observations)),
                  covariate.labels = c("Treatment: Merit", "Treatment: Patronage", "Intercept (Reference: Random)"),
                  
                  #star.cutoffs = c(0.05, 0.01),
                  #float.env = 'sidewaystable',
                  keep.stat = c("n"),
                  notes = NULL,
                  notes.align = 'l')

write_latex(table[-c(10, 11, 12, 18, 21, 24, 30)], note_text, './outputs/tables/table_a12.tex', .8)


